Optimal and hysteretic fluxes in alloy solidification: Variational principles and chimney spacing 
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We take a numerical approach to analyze the mechanisms controlling the spacing of chimneys - channels 
devoid of solid - in two-dimensional mushy layers formed by solidifying a binary alloy. Chimneys are the 
principal conduits through which buoyancy effects transport material out of the mushy layer and into the liquid 
from which it formed. Experiments show a coarsening of chimney spacing and we pursue the hypothesis that 
this observation is a consequence of a variational principle: the chimney spacing adjusts to optimize material 
transport and hence maximize the rate of removal of potential energy stored in the mushy layer. The optimal 
solute flux increases approximately linearly with the mushy layer Rayleigh number. However, for spacings 
below a critical value the chimneys collapse and solute fluxes cease, revealing a hysteresis between chimney 
convection and no flow. 

PACS numbers: 47.20.Bp, 47.20.Hw, 05.70.Ln, 47.54.-r 
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Variational principles constitute a cornerstone of physics 
because the trajectory of a system is determined from the ex- 
tremum of the action. A common example in classical physics 
is an action defined as the time integral of the Lagrangian. 
However, variational principles for nonlinear dissipative sys- 
tems constitute a topic of long standing debate because the 
non-conservation of phase space volume implies such systems 
are not Hamiltonian [e.g.Q]]. Successful examples include tur- 
bulent Rayleigh-Benard convection, where a variational ap- 
proach yields bounds on the heat flux that compare favorably 
with scaling arguments J2l; a similar approach has been ap- 
plied to shear driven turbulence In this Letter we consider 
how a variational principle can be applied to describe convec- 
tion in a mushy layer, a reactive porous medium formed dur- 
ing solidification of a binary alloy |4|]. In addition to shedding 
light on the dynamics of nonlinear dissipative systems, this 
problem has direct applications in geophysical, geological and 
industrial settings. For example, mushy-layer convection is 
principally responsible for brine drainage from young sea ice 
and the consequent buoyancy forcing of the polar oceans @]. 

Under common growth conditions morphological instabil- 
ity of the solid-liquid interface generates a mushy layer: a 
reactive porous medium of solid dendrites bathed in con- 
centrated fluid. The interstitial fluid can become convec- 
tively unstable resulting in buoyancy-driven convection within 
the mushy layer |6fl. Convection drives flow of solute de- 
pleted/enriched fluid into regions of high/low solute concen- 
tration, leading to local dissolution/growth of the solid matrix 
because fluid in the interstices adjusts to maintain local ther- 
modynamic equilibrium. Regions of low solid fraction have 
high permeability, and hence flow focussing accelerates the 
growth of the instability. The nonlinear growth of this insta- 
bility leads to the formation of channels of zero solid fraction, 
or chimneys, which form the principal conduits for drainage of 
solute from the layer. Experiments show that under a constant 
solidification rate, the chimneys are regularly spaced [e.g.0], 
whereas during growth from a fixed temperature surface the 
mean spacing between chimneys increases over time as the 
mushy layer thickens . 



The onset of convection and local dissolution is predicted 
by linear and weakly nonlinear stability analyses (reviewed 
in ylla] ), but after chimneys form a different theoretical ap- 
proach is required to account for the combination of porous 
medium flow in the mushy region and pure liquid flow in the 
chimney. Previous analyses treated either an isolated chim- 
ney |@] or modelled dynamics that arise with a periodic array 
of chimneys of imposed spacing ifiol [Till . These require the 
areal number density of chimneys to be specified a-priori, and 
thus any subsequent prediction of solute fluxes relies on an in- 
dependent theoretical prediction of the spacing of chimneys. 

We hypothesize that the chimney spacing adjusts to opti- 
mize drainage of potential energy from the mushy layer, and 
thus the system dynamics are determined by a variational prin- 
ciple that yields optimal solute fluxes. The resulting properties 
are determined numerically for two-dimensional steady-state 
solidification and we use this to reconcile behavior observed 
during transient growth. Fig.[T]describes the two-component 
mixture of liquid concentration C and temperature T that is 
translated at a velocity V between hot and cold heat exchang- 
ers. Solidification depletes the liquid of solute, reducing the 
density and providing the buoyancy that drives convection. 
We investigate the behavior of a periodic array of chimneys 
within this system. For a given chimney spacing I, we calcu- 
late the resulting solute flux. 

We employ so-called "ideal" mushy layer theory which de- 
scribes conservation of heat and solute, along with incom- 
pressible Darcy flow. We draw upon an analysis of boun dary 
conditions and their implications developed previously (Kj- 
[l3ll . and combine this with a fully time-dependent treatment 
and the hypothesized variational principle to reveal the new 
results presented here. Within the mushy layer T and C are 
coupled by local thermodynamic equilibrium and hence lie on 
the liquidus curve T = T L (C) = T E + T(C - C E ), where 
r is constant, so that the local dimensionless temperature is 
= [T — T L (C )] /TAG = (C - C ) /AC where C is 
the concentration in the liquid layer and AC = Co — Ce. 
We solve for the dimensionless temperature 8 and solid frac- 
tion <t> and calculate the Darcy velocity u = (tp x , —ip z ) by 
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z = h(x, t) 



FIG. 1. (Color online) A two dimensional mushy layer formed be- 
tween two heat exchangers pulled at velocity V aligned with the 
gravitational acceleration gk. The overlying liquid has constant far- 
field temperature and concentration Co, and the mush solidifies 
at eutectic temperature Te and concentration C'e at the lower bound- 
ary z = 0. We assume a periodic array of chimneys, of dimensional 
width 2a(z, t) and imposed spacing I, and exploit symmetry to solve 
for the properties of a mushy layer of thickness h(x, t) within the 
dashed outline. The specific heat capacity c p and thermal diffusivity 
K are assumed constant across solid and liquid phases, and the fluid 
has dynamic viscosity p, and density pog/3(C — C'e) for constant ha- 
line coefficient /3 and reference density pa. This occurs in aqueous 
NH4CI solidified from below and the dynamics and thermodynamics 
are ostensibly the same as ice forming above salt water. 



ditions at the right hand boundary of the domain x = A. 

The boundary conditions at the chimney wall x = a(z,t) 
play a key role in describing the flow. Chimneys are narrow 
(a <C 1) so they can be represented by singular interface con- 
ditions at x = 0. Lubrication theory applied to the flow in the 
narrow chimney yields the mass flux condition 
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where the pre-factor for the forcing has been calculated from 
a quadratic Polhausen approximation iflolfTlll . Balancing the 
heat flux conducted into the chimney with that advected along 
the chimney yields 
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The chimney wall is a free boundary with net outflow and 
radius a(z, t) determined from the condition lfl3ll 
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generating a vorticity equation for the dimensionless stream- 
function ip, assuming that the fluid density depends linearly 
on concentration and that the mushy layer has permeability 
TI = n (l — 0) 3 . Velocities, lengths and times are scaled by 
V, It = k/V and tx = k/V 2 respectively from which we 
obtain six dimensionless parameters governing the system, 
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The mushy layer Rayleigh number R m describes the ratio of 
buoyancy to dissipation, and the Darcy number Da charac- 
terizes the mushy layer permeability. The Stefan number S, 
concentration ratio C and scaled temperature 9^ characterize 
the imposed thermodynamic conditions. 

Rather than solving directly for the overlying fluid layer, 
we apply a boundary layer approximation to describe its in- 
fluence on the mushy layer. We assume constant pressure at 
the mush-liquid interface (12], and that in the absence of so- 
lutal diffusion the fluid region has uniform concentration Co 
away from plumes exiting the mushy layer fioll . The position 
of the mush-liquid interface is determined by the condition of 
marginal equilibrium 9 = at z = h, and hence continuity of 
salinity and normal heat fluxes give = Oand n ■ VT|+ = 
at z = h. Applying a boundary layer approximation that bal- 
ances advection and diffusion of heat across isotherms of cur- 
vature V • n yields 



n • V0 = 0oo [V • n - (u - k) • n] . 



(2) 



The lower boundary z 
eutectic temperature (0 



is impermeable and fixed at the 
— 1), and we apply symmetry con- 



The system, including boundary conditions ©-(O, was 
integrated numerically using second-order finite differences, 
with heat and concentration equations treated using semi- 
implicit Crank-Nicolson time-stepping. Elliptic equations for 
xjj and were solved using multigrid iteration lfl4l IT5I1 . Fi- 
nally, the chimney radius and mush-liquid interface position 
h(x, t) were updated using relaxation. The chimney radius 
was treated as a free boundary and updated at each spatial 
grid-point to reduce the error in (O. The boundary layer ap- 
proximation d2J leads to an unstable scheme for the corre- 
sponding free boundary problem for h(x, t). Hence, we en- 
force a one parameter shape 

h = hi — t/j c [1 — cosh /i (A — x)} I (/isinh /iA) , (6) 

where tp c is the streamfunction value at x = (0, h) [e.g.. [Toll . 
To remove a temperature singularity at the chimney top this 

— 2/3 

shape has a thermal boundary layer of width = R m 
The parameter hi is adjusted to minimize the residual in sat- 
isfying (|2]) in a least squares sense, with = enforced at 
z = h. Importantly, the time-dependent initial value prob- 
lem was integrated to a steady state, for imposed values of 
the chimney half-spacing A = IV/2k . The initial conditions 
were given either by a similarity solution with no fluid flow ^ 
or by continuation from a previous steady state solution with 
different parameters. An arc length continuation scheme was 
also used to provide an alternative confirmation of the steady- 
state solution branches JTF 



We investigate the influence of chimney spacing and con- 
vective strength on the mushy layer dynamics, using solutions 
for a range of A and R m with C, S, 9 X and Da held fixed. First 
consider the variation of the chimney spacing wavelength A 
with the Rayleigh number held fixed at R m = 40, noting 
that qualitatively similar behavior is observed for other val- 
ues 27 < R m < 60. The solute flux F$ from the mushy 



3 




X = IV/2k X = IV/2k 



FIG. 2. (Color online) (a) Steady-state solute fluxes Fs (A) exiting a 
mushy layer vary with the chimney spacing A. Two branches of solu- 
tions are determined and exhibit hysteresis; the upper branch results 
from convection with chimneys (red crosses) and the lower branch 
results from a state of no flow and hence -F!s(A) = (blue squares). 
The upper branch exists for A > A s and the lower branch is unstable 
for A > An with red and blue dashed lines indicating the solution tra- 
jectories at these points. A maximal solute flux Fso is attained for a 
chimney spacing A = Ao, with weaker solute drainage at A > Ae>. 
The inset shows detail of the upper stable branch (red curve) and in- 
termediate unstable branch (black dashed curve) in the vicinity of 
the stabilization point confirmed by arc length continuation. The 
calculations are for R m = 40, C = 15, «S = 5, #oo = 0.4 and 
Da = 5 x 10~ 3 for consistency with previous studies ITTIl which use 
properties for aqueous NH4CI. The optimal solute flux is approxi- 
mated by ( 1718b . with 7 = 0.03 and R c = 20 for these parameters, 
(b) Stability curves tracing the variation of A s (red points) and \ u 
(blue points) with R m , with all other parameters held fixed. In I only 
the lower branch is stable yielding no flow, and in II only the up- 
per branch is stable yielding convection with chimneys. Hysteresis 
is observed with two steady states in III. Both no flow and chimney 
convection states are unstable in IV, and we observe a state of weak 
convection with no chimney. 



layer is shown as a function of chimney spacing A in fig.[2a). 
There are two steady state branches, a lower branch corre- 
sponding to a state of no flow, and an upper branch describ- 
ing convection with chimneys. Depending on the choice of 
initial conditions, hysteresis is found with one of two stable 
steady solutions over a range A s < A < X u . A state of no 
flow remains stable for A < X u , but becomes unstable for 
A > X u with the solution evolving in time to the upper branch 
of chimney convection. If we start on the upper branch and 
reduce A then chimney convection remains stable for A > X s , 
but when A < A s chimneys collapse, returning the system to 
a state of no flow. Fig. [2|b) traces the stability boundaries 
X s and A u versus R m , and identifies regions of phase space 
with no flow (I), chimney convection (II) and both steady 
states (III). Hence, starting in a state of chimney convection 
in region III and reducing A the system crosses the stability 
boundary A s (R m ), the flow is stabilized and chimney convec- 
tion ceases. For R m < 27, the stability curves cross and the 
nature of the solution changes. Additional calculations indi- 
cate a state of weak convection with no chimneys, and hence 
Fs=Q, observed in region IV, where both chimney convection 
and no flow states are unstable. 

Because there are always sufficiently large wavelengths A 
available to trigger the instability of a state of no flow, we ex- 
amine the the upper solution branch to find that Fg(X) has a 




0.2 0.4 0.6 0.8 




x 




FIG. 3. (Color online) Comparison of mushy region profiles at a 
long-wavelength chimney spacing A = 1.0 > Ao (a,b), with those at 
the optimal chimney spacing A = Ao = 0.46 (c,d). The temperature 
9 is indicated by the color scale in (a,c) with isotherms shown as solid 
black curves. Solid fraction <f> is shown by the color scale in (b,d), 
constant <j> contours as black curves, and Darcy velocity streamlines 
as magenta curves. Other parameters are identical to fig. [2] 



maximum at some critical wavelength X — Xo- Hence, an op- 
timal solute flux can be attained by varying the chimney spac- 
ing (Fig. [2^). The solute flux weakens at large wavelengths 
A > Ao, which can be understood by considering examples 
of the mushy layer properties at different chimney spacings. 
Figs. [3ja,fe) show profiles of steady state mushy layer tem- 
perature, solid fraction and streamlines of Darcy velocity for 
A = 1.0 at A > Xq- At this large wavelength, approximately 
half of the mushy region is well drained by streamlines en- 
tering at the upper boundary and exiting through the chimney 
at x = 0. However, there is a large nearly stagnant region 
away from the chimney, suggesting an explanation for the ob- 
served inefficient drainage for large chimney spacings. Com- 
pare this to the corresponding profiles for the optimal config- 
uration at A = Ao (Figs. |3J:,d) where the streamlines show 
efficient drainage via convective cells of order one aspect ra- 
tio. Thus, rather than drainage rate being controlled by buoy- 
ancy driven flow in the chimney, the optimal solute flux is 
controlled by the efficiency of convection within the mushy 
region. The temperature and solid fraction have qualitatively 
similar structure for both wavelengths, with significant hori- 
zontal variation of the latter leading to inhomogeneity in the 
concentration of the final material. 

Having determined the system properties for a range of A, 
we now apply the variational principle to select a preferred 
value of chimney spacing with maximal solute flux at A = Ao 
and calculate how the system varies with R,„. The R m ^> 1 
simulations show that the optimal solute flux Fso increases 
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approximately linearly with R m suggesting the approximate 
scaling laws 

F s = Rm < R c (7) 

F s ~ 7 (R m -R c ) R m > R c , (8) 

for some constants 7 and R c that depend on the other parame- 
ters imposed on the system. For R m ^> 1 Eq. (O implies that 
the dimensional solute flux 

F s ~ JPogPiCo - C E ) 2 H /fi, Rm » 1, (9) 

is independent of both the thermal diffusivity k and the so- 
lidification rate V, This is consistent with the rate of solute 
transport being controlled by the large scale convective flow, 
independent of any effective transport induced by molecular 
diffusion. As a point of comparison, the heat flux in turbulent 
Rayleigh-Benard convection is also predicted to be asymptot- 
ically independent of k in Kraichnan's ultimate strongly con- 
vective regime lfl7tl . 

The optimal chimney spacing Xo and resulting mushy layer 
depth ho both decrease as the Rayleigh number increases but 
the aspect ratio Xo/ho asymptotes to a constant value for 
R m 3> 1 . The stronger flow at larger R,„ generates a thinner 
mushy layer, but the most efficient solute drainage is given 
by order one aspect ratio convective cells. This behavior is 
consistent with the constant mean aspect ratio observed in the 
transient phase of enthalpy method simulations flUl . 

These results embolden us to suggest explanations for phe- 
nomena observed during transient growth, such as growth 
from a fixed temperature boundary. Experiments show that, 
as the mushy layer thickens over time, extinction of convec- 
tive flow in some of the chimneys leads to an increase of the 
mean spacing of chimneys [511 . This coarsening may be con- 
sistent with the dynamics of optimal chimney spacing which 
we find has constant aspect ratio Xo/ho for R m 3> 1. This 
is consistent with the mean spacing of chimneys increasing 
with mushy layer depth during transient growth. Moreover, 
during transient growth, the extinction of flow in certain con- 
vective channels may be consistent with the flow stabilization 



found here for A *C h. Because the mean depth h increases 
with fixed chimney spacing A during transient growth, the 
aspect ratio X/h decreases until it triggers a stabilization of 
convection and extinguishes flow in a selection of the chim- 
neys. Taken together this offers a possible explanation for the 
observed mechanisms of the coarsening of chimney spacing 
as h increases. Comparison with previous work 151 1131 sug- 
gests that the scaling (O may also be of relevance for transient 
growth at small concentration ratios (C <C 1). In particular, 
consistent with experiments in a finite geometry Jsj], the solute 
flux (O predicts that the concentration of the liquid region will 
change approximately linearly in time. This would provide a 
simple parameterization of brine drainage from growing sea 
ice for use in large scale models without having to resolve 
natural horizontal variations in sea ice structure. 

In summary, we have numerically analyzed strongly nonlin- 
ear convection in a solidifying mushy layer with a periodic ar- 
ray of chimneys with spacing A. By varying A, we have shown 
the existence of an optimal chimney spacing Xo that maxi- 
mizes the solute flux from the mushy layer and hence also the 
rate of removal of its potential energy. This Ae> yields convec- 
tive cells of order one aspect ratio thereby efficiently draining 
the mushy layer, with weak flow for A ^> Xo- For A <C Xo 
there is stabilization, so that chimney convection cannot be 
supported for spacings smaller than a Rayleigh number de- 
pendent critical value, which suggests a method to suppress 
chimney formation in engineering applications. Steady states 
of chimney convection and no flow show hysteretic behavior. 
These mechanisms are consistent with dynamics controlled 
by a variational principle, with the spacing of chimneys ad- 
justing to optimize the rate of release of potential energy from 
the mushy layer, and facilitate the most efficient route towards 
thermodynamic equilibrium. 
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